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ABSTRACT 

The back-reaction of baryons on the dark matter halo density profile is of great in- 
, terest, not least because it is an important systematic uncertainty when attempting 

Qh' to detect the dark matter. Here, we draw on a large suite of high resolution cos- 

Q I mological hydrodynamical simulations, to systematically investigate this process and 

5^ . its dependence on the baryonic physics associated with galaxy formation. The effects 

• of baryons on the dark matter distribution are typically not well described by adi- 

I abatic contraction models. In the inner ten per cent of the virial radius the models 

are only successful if we allow their parameters to vary with baryonic physics, halo 
mass and redshift, thereby removing all predictive power. On larger scales the pro- 
files from dark matter only simulations consistently provide better fits than adiabatic 
contraction models, even when we allow the parameters of the latter models to vary. 
I The inclusion of baryons results in significantly more concentrated density profiles if 

, radiative cooling is efficient and feedback is weak. The dark matter halo concentration 

^T) • can in that case increase by as much as 30 (10) per cent on galaxy (cluster) scales. 

The most significant effects occur in galaxies at high redshift, where there is a strong 
anti-correlation between the baryon fraction in the halo centre and the inner slope of 
' both the total and the dark matter density profiles. If feedback is weak, isothermal 

inner profiles form, in agreement with observations of massive, early-type galaxies. 
However, we find that AGN feedback, or extremely efficient feedback from massive 
stars, is necessary to match observed stellar fractions in groups and clusters, as well 
K> , as to keep the maximum circular velocity similar to the virial velocity as observed 

}_j ■ for disk galaxies. These strong feedback models reduce the baryon fraction in galaxies 

I by a factor of 3 relative to the case with no feedback. The AGN is even capable of 

reducing the baryon fraction by a factor of 2 in the inner region of group and clus- 
ter haloes. This in turn results in inner density profiles which are typically shallower 
than isothermal and the halo concentrations tend to be lower than in the absence 
of baryons. We therefore conclude that the disagreement between the concentrations 
inferred from observations of groups of galaxies and predictions from simulations that 
was identified by Duffy et al. (2008) is not alleviated by the inclusion of baryons. 

Key words: galaxies: clusters: general - galaxies: haloes - dark matter - cosmology: 
theory methods: A'^-body simulations 
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1 INTRODUCTION 

Increasingly powerful computers, combined with ever more 
efficient Af-body codes, have permitted accurate numer- 
ical tests of structure f ormation in a Co l d Dark Mat- 
ter (CDM) universe (e.g. iKlvpin et allll996l : iNavarro et~al] 



19971: 



20051 : 



Moore et al]|l999l: iBuUock et al.l 



Diemand et al.l l2008l : iGao et al.l 



20011 : ISpringel et al.l 
20081 '). It is now an 



established prediction that CDM haloes form a cuspy mass 
distribution t hat is close to a 'universal profile', indepe ndent 
of halo mass (jNavarro et al.|[l996l : iNavarro et al.|[T997l . here- 
after NFW). In detail, however, the haloes are not strictly 
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self-similar (|Navarro et al.|[201Cll V Foremost amongst these 

results has been the existence of power law relationships in 

the ph ase-space density of halo es across several decades in 

radius (iTavlor fc Navarrcll200ll l and the halo concentration- 

\ " 1 

mass relation over five decades in mass (Bullock ct al. 2001^; 
iNeto et al, 2007: Maccio et al,.200a ; .Duffv ct al, 2008 ). High 
resolution iV-body simulations have also shown that DM 
haloes are not smooth, but contain self-bound s ubstructures 
that have survived the a ccretion process (e.g. iMoore et all 
ll999l : lKlvpin et al.lll999h . More recently, it has been shown 
that substructures themselves contain s ubstructures, a pat- 
tern that shows no signs of abating jStadel et al.l |2009| : 
ISprineel et aDl2008l l. 

These power laws and embedded substructures are the 
product of an approximately scale-free system that inter- 
acts only via gravity and, for scales large enough that grav- 
ity is the only dynamically significant force, can be viewed 
as an accurate approximation to structure formation. On 
smaller scales, other processes significantly modify the struc- 
ture of the baryons, which can subsequently affect the DM 
in the densest regions. In particular, it has long been known 
that a galaxy forms because of the ability of its pro genitor 
gas t o cool efficiently via emission of radiation (e.g. IHov13 
1 19531 ). Such a process introduces a characteristic physical 
scale into t he system and thus breaks the self-similarity of 
structures (| White fc Reeslll978h . The net effect is that the 
cooled gas and stars then cluster on smaller scales than the 
DM, thereby deepening the potential well of the system and 
causing the DM to orbit closer to the centre. This will modify 
a number of A'^-body simulation predictions on small scales, 
in a way that depends on the efficiency of galaxy formation. 

This modification of the DM distribution is usually 
parametrised in the form of adiabatic invariants, hence the 
term 'adiabatic contraction' that is used to describ e the 
overall effect, first co nsidered by lEggen et al.l (| 19621 ) and 
IZeldovich et all ljl980l l. The adiabatic parameter considered 
in the first gener ation of adiabatic contr action models for 
collapsing haloes (|Blumenthal et al]|l986l ) was the product 
of radius with the mass internal to this radius, M(< r)r. 
iGnedin et al.l (|2004l 'l. henceforth G04, found that such a sim- 
ple scheme, only strictly applicable for particles on circu- 
lar orbits, did not provide a good description of what was 
happening in their hydrodynamic simulations of groups and 
clusters at 2: « and galaxies at 2: > 3. Better agreement 
was found when they modified the parameter to use the 
orbit-averaged radius, f, that is a power-law function of ra- 
dius. The work of G04 was extende d to galaxy scales at low 
redshift: bv lGustafsson et all (j2006l ). who modified the slope 
and normalisation of this power law distribution to provide a 
better description of their results at this lower mass range. 
The compression of the DM by the infalling baryons has 
also been extensively studied by purely analytic methods 

iKlar 



Sellwood fc McGaughl [2OO5I : [ Cardone fc Serend [2OO5I : 



Miic"ketll2008l '). 

In addition to the contraction of the DM halo, the con- 
centration of baryon ic mass at its centre also leads to a more 
spherical structure jKazantzidis et al]|2004l 'l. This effect is 
due to the more spherical gravitational potential formed by 
the baryons and it is this that induces an adiabatic decrease 
in the eccentricity of the orbits of the DM and stars, rather 
than through the destruction of box orbits jDebattista et all 
|2008| ). Additionally, in the presence of a gaseous disc, in- 



falling subhaloes with low orbital inclinations with respect 
to the disc will experience significant disruption due to the 
high local density of the baryons. This will lead to the accu- 
mulation of stars, and DM, in the same plane, f orming both 
a thin stellar and a DM disc (|Read et al.ll2008l ). 

With the addition of baryons to the simulation, one 
can expect further effects due to the tr ansfer of angular 
momentum from infalling sate llites (e.g. iDebattista et al] 
l2008l : iRomano-Diaz et al.ll200§ '). However, the infiuence of 
baryons on these smaller, less well-resolved objects, are by 
no means certain. For example, the substructure may ex- 
hibit greater resistance to tidal disruption due to its deeper 
potential well in the presence of baryons, while at the same 
time the substructure will typically su ffer from increase d 
dynamical fricti on with the main h alo (| Jiang et al.ll200^ '). 
Recent work by iMaccio et al.l (|2006l ) found a factor of two 
increase in the number of surviving satellites in a galaxy- 
sized halo, in a hydrodynamical simulation relative to a DM- 
only si mulation. This increas ed survival rate in the inner 
regions (|Weinberg et al.ll2008l ) may enable satellite infall to 
partly counteract the contraction of the DM halo, as these 
objects can efficiently transfer angular momentum to the 
inner parts of the DM halo, instead of being tidally dis- 
ru pted at larger radii ( f or recent considerat ions of this effect 



IPedrosa et al.ll2009l . IPedrosa et allboiol and lAbadi et all 
120091 ). The erasure of the central DM cusp in galaxies b y in- 



falling satellites was prop osed by El-Zant et al] (|200ll ) and 
for the case of clusters bv lEl-Zant et al. (2004). 

The effect of the baryons on the inner DM density pro- 
file is of particular interest, as it is currently a significant 
source of uncertainty when making predictions for DM de- 
tection experiments. For example, the expected 7-ray sig- 
nal f rom self-annihi lation of potential DM candidate parti- 
cles jSpringel et aLlfeoOS. ). which may become detectable in 
the near future using 7-ray observatories such as fermQ, is 
proportional to Pdm- We further note that the total mass 
density profile is also of interest, as it is this quantity that 
is determined from strong lensing studies and is relevant for 
direct comparison with such observations (for a review of 
the m ethodologies and uses of strong lensing, see iKochanekl 
120061 ). 

The incorporation of baryons in a simulation is a chal- 
lenging theoretical undertaking as the relevant scales of the 
physical processes are rarely resolved in a cosmological sim- 
ulation and approximations are therefore necessary. The 
method by which stellar feedback, for example, is incorpo- 
rated can have a large infiuen ce on the resulting baryonic dis - 
tribution of the galaxy (e.g. iDalla Vecchia fc Schavell200"8l ). 
For this reason, we have attempted to probe these effects in 
a systematic way by examining the DM halo density profiles 
from galaxies, groups and clusters, drawn from an extensive 
series of cosmological simulations with that incorporate a 
wide variety of prescriptions for baryonic processes. These 
form a subset of the larger, overall simulation series known 
as the Overwhelm ingly Large Simulations project (OWLS; 
ISchave et al]|2010l ). 

Our work has a number of novel features that allow us 
to extend recent studies that have probed the dynamical re- 
sponse of the DM halo profile to the presence of baryons. For 
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example, while iGustafsson et all (|2006l ) examined a small 
number of galaxy-sized haloes at z — 0, we provide results 
on both group and cluster mass scales, with significantly 
better statistics. G04 examined the DM profiles in a simi- 
lar mass range to our own, albeit limited to higher redshift 
for galaxy scales, and with similar mass resolution. They 
demonstrate the significant impact that metal enrichment 
and stellar feedback have on the overall gas distribution, al- 
though they did not vary the mod els systematically. The re- 
cent work of IPedrosa et al.l (|2009l l investigated the effects of 
physics on galaxy scales by simulating a single halo. Here we 
provide a more statistically-robust set of results, due to the 
large number of haloes we have resolved. (We note that a di- 
rect comparison with their z = galaxy is unfortunately not 
yet possible and must await further simulations). In sum- 
mary, the large dynamic range in mass, from dwarf galaxies 
to clusters, across a wide variety of physics implementations 
and with high statistical significance, is a significant advance 
in the study of baryons in DM haloes. 

The rest of the paper is organised as follows. In Sec- 
tion[2]we describe the range of simulations used in this study, 
along with the differences between implementations of gas 
physics. In Section [3] we investigate the distribution of the 
baryons within these haloes as a function of the physics im- 
plementation. In Section |4] we compare our haloes with both 
strong gravitational lensing observations at high redshift and 
stellar mass fractions at the present day. We discuss the 
implications of our findings for structure formation in the 
CDM model, with the proviso that certain simulations can 
fulfil some observations but none can be used at all mass 
scales. The analysis of the haloes formed in these simula- 
tions, via their density profiles, is discussed in Section[5l We 
then attempt to parametrise the DM density profile with 
existing theoretical profiles, and the total matter profile by 
comparing non-parametric measures of halo concentrations, 
in Section [6] In Section [7] we test whether we can mimic 
the infiuence of the b aryons by making use of the adiabatic 
contraction models of lBlumenthal et al.l (| 19861 ) and G04. Fi- 
nally, we summarise our work in Section [S] 



2 SIMULATIONS 

This work is based on OWLS; a series of high-resolution 
simulations of cosmologic al volumes with di ffering subgrid 
physics implementations (|Schave et al.|[201ol ). These simu- 
lations were run using a modified version of GADGET-3, it- 
self a modified vers ion of the publicly- available GADGET-2 
code (|Springelll2005h . GADGET calculat es gravitational forces 
using the Particle-Mesh algorithm fe.g. lKlvpin fc Shandariiil 



1 19831 ) on lar ge scales, supplemen ted by the hierarchical 
tree method ( Barnes fc Hull [19861 ) on smaller scales. Hy- 



drodynamic forces are computed us ing the Smooth Particle 
Hydr odyna mics (SPH) formalism jMonaghan fc Lattanzid 
Il98g : Sprin gel fc Hernauist|[200^ '). 

The production simulations available in OWLS model 
cubes of comoving length 25 and 100/i^^Mpc with 512'^ 
gas and 512'' DM particles. For all runs, glass-like cosmo- 
logical initial conditions were generated a.t z = 127 using 
the Zel'dovich approximation and a transfer function gen - 
erated using cmbfast (v. 4.1: ISeliak fc Zaldarriaga|[l996l ). 
Cosmological parameters were set to the best-fit values 



from the 3rd year Wilkin son Microwave Anisotropy Probe 
data ISpergel et all 120071 ). henceforth known as WMAP3. 
Specifically, the values for [Qm, f^b, i^A, h, as, ris] were set 
to [0.238, 0.0418, 0.762, 0.73, 0.74, 0.95]. The primordial 
mass fraction of He was set to 0.248. All runs that used the 
same box size had identical initial conditions. 

We have analysed a subset of 6 OWLS runs for each 
box size: a DM-only run (hereafter labelled DMONLY) and 
5 additional runs that also follow the baryonic component 
with hydrodynamics, gas cooling and star formation. For 
DMONLY, the particle masses are m = 7.7 X 10^ and 
4.9 X 10*/i"^Mo for the 25 and WOh'^Mpc runs, respec- 
tively. For the runs with baryons, each particle was split in 
two in the initial conditions, with the mass shared according 
to the universal baryon fraction, /b"'^ = fib/^^m = 0.176, 
such that the DM (gas) particle mass is 6.3 (1.4) X 10® and 
4.1 (0.86) X 10* h-'^MQ for the 25 and 100 h'^Mpc runs, re- 
spectively. The 25/i~^Mpc simulations were run to z = 2 
while the 100/i~^Mpc simulations were run all the way to 
2 = (our results will focus on data from these redshifts). 
At early times the softening length was held constant in co- 
moving co-ordinates at 1/25 of the initial mean inter-particle 
spacing. Below z — 2.91 the softening length was held fixed 
in proper units (thus, the 2 = values were 0.5 /i~^kpc and 
2/i~^kpc respectively). 

2.1 Baryonic physics 

For our purposes two main aspects of the baryonic physics 
have been varied: radiative cooling (with and without metal 
lines) and feedback (from supernovae and accreting super- 
massive black holes). The modification of the baryon distri- 
bution will be primarily determined by these two competing 
effects. Within OWLS other physics prescriptions have been 
varied, such as the choice of stellar initial mass function for 
example, however these will at most have a secondary im- 
pact on the baryon distribution. We use 'PrimC' and 'ZC' 
to label two different approaches to cooling, denoting the 
absence- or inclusion of metals to the cooling rate respec- 
tively. For the supernova feedback, we use 'NFB', 'WFB' 
and 'SFB' to identify whether stellar feedback is absent, 
weak or strong. In addition, we have analysed a simulation 
that includes feedback due to accreting supermassive black 
holes, associated with Active Galactic Nuclei; this run is la- 
belled 'AGN'. Table [T] lists the runs wi th baryonic physic s, 
along with their original OWLS names jSchave et al.ll201Gl ). 
and highlights their key differences. We only provide brief 
descriptions here, for further details please see lSchave et al.l 
(|2010l ) and the references given below. 

In all the simulations with baryons, radiative cooling 
(and heating) rates are implemented element-by-element us- 
ing t ables generated wit h the CLOUDY radiative transfer 
co de jFerland et al.lll998l ). following the method described 
m |w lersma et all (|2009l ) We assume coUisional equilibrium 
before reionisation (z > 9) and photoionisaton equilibrium 
afterwar ds, in the presence of a n evolving UV/X-ray back- 
ground (|Haardt fc Madaull200ll ). In the 'PrimC' runs cool- 
ing rates were computed assuming primordial abundances. 
In those labelled 'ZC we also track line emission from nine 
metals: C, N, O, Ne, Mg, Si, S, Ca and Fe. Both types of 
simulation also include cooling via free-free Bremsstrahlung 
emission and Compton cooling due to interactions between 
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Table 1. A list of all simulations, the corresponding names in the OWLS project jSchave et al.|[2"oi(J) and the differences in the included 
subgrid physics. 



Simulation 


OWLS name 


Cooling 


Supernova 


AGN 








Feedback 


Feedback 


PrimCJ^FB 


NOSN.NOZCOOL 


Primordial 


None 


None 


PrimC.WFB 


NOZCOOL 


Primordial 


Weak 


None 


ZC.WFB 


REF 


Metal 


Weak 


None 


ZC.SFB 


WDENS 


Metal 


Strong 


None 


ZC.WFB.AGN 


AGN 


Metal 


Weak 


Yes 



the gas and the cosmic microwave background. For full de- 
tail s of the implementatio n of gas cooling in the simulations, 
see lWiersma et al.l (|2009l ). 

When gas is sufRciently dense it is e xpected to b e 
multiphase and unstable to star formation l|Schave|[2003 ). 
Our simulations lack both the resolution and the physics 
to model the cold interstellar gas phase and we therefore 
impose an effective equation of state, P oc p^°'^, for den- 
sities that exceed our star formation threshold of riH > 
0.1 cm~^. The normalisation is fixed such that P/k — 1.08 x 
10^ K cm at the star formation density threshold. The in- 
dex is set to 7efi =4/3, which has the advantage that both 
the Jeans mass and the ratio of the Jeans length to the SPH 
smoothing length are indepe ndent of density, thus suppress - 
ing spurious fragmentat i on JSchave fc Dalla Vecchial l2008h . 
ISchave fc Dalla Vecchial (|2008l ) showed how the observed 
Kennicutt- Schmidt law, E* = ^(^g/l Mopc"^)" can be 
analytically converted into a pressure law, which can be 
implemented directly into the simulations. This has the 
advantage that the parameters are observables and that 
the simulations reproduce the input star formation law ir- 
respective of the assumed equation of state. We use the 
ISchave fc Dalla Vecc hia (,200& ) m ethod, setting A = 1.515 x 
10"* Mq yr"^ kpc"^ and n = 1.4 (jKennicuttll 19981 . note that 
we renormalised the observed relation for a Chabrier IMF). 

These star particles are assumed to be sim ple stellar 

J lobulations with an initial mass function given bv lChabrieJ 
2003). The metal abundances of the stellar particles are in- 
herited from their parent gas particles. The stellar evolution, 
and the subsequent release of metals from the star particle, 
depends on this metallicity (|Portinari et al.l 1 19981 ; iMarigol 
l200ll : lThielemann et aLll2003l ). The delayed release of 11 in- 
dividual elements (the ones used for the calculation of the 
cooling rates) by massive stars, from Type la and Type II 
supernovae, as well as AGB stars, is tracked. Every time 
step they are shared amongst neighbouring gas particles ac- 
cording to the SPH interpolation scheme. For further details 
please see I Wiersma et all llooi). 

Stars formed in the simulations labelled 'WFB' and 
'SFB' inject energy into local gas particles. This energy is in 
kinetic form, i.e. nearby gas particles are kicked away from 
the stars. On average, each newly formed star particle kicks 
rj times its own mass, where 77 is the dimensionless mass 
loading parameter, by adding a randomly oriented velocity 
i)w to the velocity of each kicked particle. The simulations 
labelled 'WFB' use rj = 2 and Uw ~ 600 kms~^, which cor- 
responds to forty per cent of the SN energy for our initial 
mass function. Further details of the feedback prescription 
can be found in iDalla Vecchia fc Schav3 l)2008h . Note that 



although we have labelled the supernova feedback in this 
simulation "weak", it is in fact strong compared to many 
prescriptions in the literature. 

The runs labelled 'SFB' implement wind velocity and 
mass-loading parameters that depend on the local density 
of the gas, — 600 kms~^ (nn/O.l cm~^)^^^ and r) — 

2 (uw/600 kms^^) ^, such that in the densest regions winds 
are launched with the largest speeds and smallest mass- 
loadings. For gas that follows the imposed effective equa- 
tion of state, this implies a wind speed that is proportional 
to the effective local sound speed, i;„ oc Cs,cff oc (P/pf^'^. 
The result is that winds in the 'SFB' model remove gas 
from higher mass systems more efficiently than in the 'WFB' 
model. Note that the WFB and SFB models use the same 
amount of SN energy. The difference in the efhciency thus re- 
sults purely from energy is distributed between mass-loading 
and wind velocity. 

Finally, in the simulation labelled 'AGN' supermassive 
black holes (BHs) are grown, with subsequent feedback, at 
the cent res of all massive haloes according to the method- 
ology of iBooth fc Schw3 (Hooi), w hich is a substantiall y 
modified version of that introduced bv lSpringel et all (|2005h . 
Seed BHs of mass msoed = 9 x 10* Mq (i.e. 10"'^ mg in 
the 100 Mpc box, where mg is the initial mass of the 
gas particles), are placed into every DM halo more massive 
than 4 x lO" Mq (i.e. 10^ DM particles in the 100 h'^ Mpc 
box). Haloes are identified by regularly running a Friends- 
of-Friends (FOF) group finder on-the-fly during the simula- 
tion. After forming, BHs grow by two processes: accretion of 
ambient gas and mergers. Gas accretion occurs at the mini- 
mum of the Eddington rate and the Bondi-Hoyle (|l944l ) rate, 
where the latter is multiplied by (nn/lO"^ cm~^)'^ for star- 
forming gas (i.e. TiH > 10~^ cm"'') to compensate for the 
fact that our effective equation of state strongly underesti- 
mates the accretion rate if a cold gas phase is present. Gas 
accretion increases the BH masses as ttibh = (1 — er)wiaccr, 
where tr = 0.1 is the assumed radiative efficiency. We as- 
sume that 15 per cent of the radiated energy (and thus 1.5 
per cent of the accreted rest mass energy) is coupled ther - 
mally to the surrounding medium. iBooth fc Schavj (|2009l ) 
showed that this model reproduces the redshift zero cosmic 
BH density as well as the observed relations between BH 
mass and both the stellar mass and the central stellar ve- 
locity dispersion. 



2.2 Halo definitions and density profiles 

The principal quantity that is extracted from each simula- 
tion is the spherically-averaged halo density profile, dom- 
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inated in all but the central regions by the dark matter. 
Haloes were first identified using the FOF algorithm, that 
links dark matt er particles using a dimensionless linking- 
length, b = 0.2 ( Davis et aDllQSSl). We then used the S UB- 
FIND algorithm (jSprineel et al.ll200ll : [Polag et al1l2009l 'l to 
decompose the FOF group into separate, self-bound sub- 
structures, including the main halo itself. Finally, a sphere 
was placed at the centre of each halo, where the centre was 
identified with the location of the minimum potential par- 
ticle in the main halo. The spheres were then grown until a 
specified mean internal density contrast was reached. 

A halo thus consists of all mass. A/a, within the radius, 
_Ra, for which 



Ma 



4 3 

-7r7?A ^ Pcrit 



(1) 



where the mean internal density is A times the critical den- 
sity, pcrit(z) = 3H{zf/8nG. If A/a = A/vir and Ra = Rvh, 
then one can compute A from the spherical top-hat col- 
lapse model. In th i s case , we adopt the fitting formula from 
iBrvan fc NormanI (|l998t ). that depends on both cosmology 
and redshift 



A = ISTr-" + 82x - S9x^ 
where x — Qin(z) — 1, 



Ho 
Hiz) 



(2) 



(3) 



and H{z) is the usual Hubble parameter for a flat universe 

H{z) = Ho y/n^{l + zy^ + QA. (4) 

In our cosmology at z = {2) A = 92.5 (168.5). We wiU 
also consider overdensities that are integer multiples of the 
critical density; A = 500 and 2500, with radii and masses 
denoted in the standard way, e.g. 7?500 and A//500 for A = 
500. 

Following iNeto et al.l (|2007l 'l. only haloes with more 
than 10'' DM particles within 7?vir are initially selected. Den- 
sity profiles are then computed using a similar method t o 
our previous work on DM-only haloes (|Duffv et al.ll200^ . 
To briefly summarise, we take 32 uniform logarithmic shells 
of width Alogj^o('') ~ —0.078, in the range —2.5 ^ 
\ogj^Q{r /Kvii) < 0. We then sum the mass within each bin 
for the three components (gas, stars and DM) and divide 
by the volume of the shell. We fit model density profiles 
by minimising the difference of the logarithmic densities be- 
tween hal o and model, ass uming equal weighting (this was 
shown bv lNeto et ahlfioOTl to be a self-consistent weighting 
scheme) . 

A key issue is the range of radii over which we can 
trust that the profiles have converged numerically. To check 
this, we have performed a convergence study, the results of 
which are presented in Appendix |X] I n summary, we show 
that the minimum radius proposed bv lPower et al.l (|2003h . 
henceforth P03, is appropriate for our simulations, even for 
the ZC_WFB model (P03 defined their resolution limit on 
the two-body dynamical relaxation timescale in the context 
of DM-only simulations). All haloes used in this work satisfy 
the criterion i?po3 < 0.05 -Rvir, so we adopt the latter value 
as the minimum radius for our density profiles. 

A selection of well-resolved haloes are singled out for in- 
vestigation of the slopes of their inner profiles (Section 15. 3p . 



For these haloes, we place even higher demands on the parti- 
cle number and ensure that haloes satisfy -Rpoa < 0.025 -Rvir 
and use the latter value as the minimum radius in this in- 
stance. In practice this criterion is applied to the DMONLY 
simulation to define a lower (virial) mass limit: at 2: = this 
is 3 X 10'^ /i"^Mq and at z = 2 it is 5 x 10^' /i"'Mq. The 
other runs are then matchecQto the DMONLY haloes that 
make this cut. All haloes used to study the inner profile have 
at least 6 x 10* DM particles within -Rvir. 

Between the two simulation sets, we have selected a 
total of 552 (204) haloes at 2 = (2) in DMONLY and of 
these, 67 (32) meet the more stringent resolution demands 
for the inner profile study. The difference in halo numbers 
in the different gas physics simulations is less than 10 per 
cent, primarily due to low-mass haloes having fewer than 
10** DM particles and hence missing the initial cut. In the 
subsequent discussion we will use the terms; Dwarf Galaxy, 
Galaxy, Group, and Gluster, to correspond to -Afvir[fe~^M0] 
in the flxed ranges [< 10'^10^^ - 10^^10'^ - 10^*, > lO'"*] 
respectively. 



3 HALO BARYON DISTRIBUTION 

The ability of baryons to cool radiatively and thereby clus- 
ter on smaller scales than the DM is a recurring theme in 
this work, leading to, potentially, numerous modifications 
of the standard ./V-body predictions. The exact distribution 
of the baryons in the halo will be sensitive to the physics 
implementation included and hence we should first quantify 
this key quantity. To that end we have plotted, in Fig. [TJ 
the baryon fraction, /b, inside the virial radius against that 
within r/Rvii — 0.05. We show the median value and the 
quartile scatter of /b, within mass bins at the different red- 
shifts. Only haloes that meet our more stringent inner profile 
criteria are shown. 

A number of interesting features are apparent from 
the figure. Firstly, we see that for Groups and Glusters 
at 2 = (bottom panel) the baryon fraction within the 
virial radius is close to the universal value in all runs except 
ZC_WFB_AGN. The other runs are approximately consis- 
tent with, though slightly hig her than, the pred iction from 
non-radiative gas simulations (jCrain et al.ll2007l ). The deep 
gravitational potential wells and long cooling times in these 
objects, lead to little gas being expelled from the halo. Only 
in ZC_WFB_AGN, which includes AGN feedback and for 
which the feedback is therefore expected to be stronger on 
group and cluster scales, does a significant amount (~ 10—40 
per cent) of gas get expelled. It should be noted, however, 
that a significant fraction of this gas will have been ejected 
from smaller systems at higher redshift. 

The inner baryon fractions, fh{r/Rvir < 0.05), typically 
decrease as the strength of the feedback increases. Note that 
ZC_WFB_AGN has an inner baryon fraction that is lower 
than the universal value. Runs with weak or no feedback 
have values that are more than twice as high in the run 



^ Haloes arc matched according to their shared DM particles; 
simulations with different physics of the same volume and reso- 
lution have identical initial conditions and particle ID lists. This 
allows us to match coincident groups of DM particles between 
runs. 
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Figure 1. The panels show the baryon fraction within the virial 
radius as a function of the baryon fraction within r < 0.05 i?viri 
at ^ = 2 (top panel) and 2 = (bottom panel). The relative 
values of these quantities gives a useful overview of the effect 
of cooling and feedback on the baryon distribution in haloes at 
different redshifts and virial masses. Each point corresponds to 
median values for haloes within a given mass bin, indicated by 
the symbol size, while the error bars indicate the quartile scatter. 
The solid vertical and horizontal lines correspond to /^"'^. The 
d ashed horizontal lin e is from the non-radiative gas simulations 
of ICrain et al. ] bOOTi ). Galaxies are more sensitive to the various 
feedback and gas cooling schemes than groups and clusters. For 
groups and clusters the sensitivity to the feedback and cooling 
is mostly limited to the inner regions but these processes affect 
galaxies all the way out to the virial radius. AGN feedback reduces 
the baryon fractions most strongly. 



with AGN feedback. The inclusion of metal-hne coohng in- 
creases the inner baryon fraction, as can be seen by compar- 
ing ZC.WFB with PrimC.WFB. 

The situation for Galaxies sd, z — 2 (top panel) is 
very different. Only the model with no supernova feedback, 
PrimC_NFB, has a median baryon fraction within the virial 
radius that is close to /b"'^; even the weak feedback model 
is capable of expelling some gas. Note that PrimC_NFB has 
a baryon fraction that is higher than predicted from non- 
radiative simulations, suggesting that gas cooling in this 
model has compensated for the heating and expansion of 
the gas due to energy transfer from the dark matter when 
the halo collapsed. In Galaxy mass haloes, the gravitational 
potential is sufHciently low that SN feedback is able to un- 
bind gas from the halo. As a result, there is a strong positive 
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Figure 2. To test the realism of our models we have compared 
the median stellar mass fractions within R500 > in units of the uni- 
versal baryon fraction, as a function of total mass within iJsoOi 
at 2 = 0. The coloured points represent the various physics pre- 
scriptions, apart fr om the single black data point at low halo 
mass which is from lConrov etal for a sample of isolated 

Lt galaxies. The vertical and horizontal error bars represent, re- 
spectively, the quartile spread and mass range within each bin. 
The solid and dashed cu r ves ill ustra te observa t ional determina- 
tions from iGiodini et al.l | |2009|) and iLin etlH ||2004|) . with the 
outer curves representing the uncertainty in the normalisation of 
their fits. It is clear that if the gas is allowed to cool by metal-line 
emission, it must be prevented from forming stars by either a su- 
pernova model that targets dense regions of the halo, ZC_SFB, or 
by feedback from AGN, ZC.WFB_AGN. We note that although 
estima tes of stellar masses are typically uncertain by a factor of 2- 
3 (e.g. iKiipcii Yoldas et"ai]|20oH iLonghetti fc SaraccdIioogI) the 
observations still favour the strong feedback prescriptions. 



correlation between fh{r/Rviv < 1) and fhir/Rvii < 0.05). 
The galaxies in both ZC.WFB and PrimCNFB have inner 
regions that are baryon-dominated (/b > 0.5). 



4 COMPARISON WITH OBSERVATIONS 

In the previous section we compared simulations run with 
a number of widely used physical prescriptions, and demon- 
strated that a large range of inner baryon fractions are po- 
tentially possible. We now determine which of these models, 
if any, are compatible with observations. We consider two 
probes that offer complementary and contrasting constraints 
on the baryon physics. First, in Section [4.11 we present a 
comparison of the predicted redshift zero stellar mass frac- 
tions to those observed in Galaxies, Groups and Clusters. 
Secondly, in Section [4. 2 1 we compare predicted total density 
profiles to those inferred for a high redshift galaxy sample 
through strong lensing measurements. 



4.1 Stellar fractions 

The stellar fraction by halo mass is well constrained obser- 
vationally therefore a comparison of median stellar fraction 
(within R500) to observation allows us to quickly determine 
the relative realism of our models. In Fig. [2] we show the 
stellar fraction, in units of M500, for haloes at z = as a 
function of mass for our full range of physics models. 
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To compare with the observation s, we hav e taken 
the group and cluste r samples of iLin et all (j2004l ) 
and iGiodini et alj (|2009h . Note that the best fitting power 
law results we take from these works do not contain the con- 
tribution from the diffuse intracluster light, and can hence 
be viewed as lower limits. This additiona l contribution likely 
ranges between 11 and 22 per cent (e.g. Zibetti et al.l 20051 : 
iKrick fc Bernstein|[2007l 'l. Following ISalogh et all ( 20081) . we 
convert the observed stellar luminosities from Lin et all 



l|2004l ) to masses, assuming the best-fit stellar ma ss-to- light 
ratio, M/Lk = 0.9. Our simulations ad opted a IChabriei 
20031') IMF, whereas iBalogh et all (|2008l ) used the iKroupa 



20011 ) IMF; we expect this difference to be unimportant, 



however, as they ar e very similar o ver th e relevant range of 
stellar masses. For iGiodini et all (|2009l ) we have adopted 
the best fitting stellar fraction for the X-ray COSMOS 
groups/poor clusters onl y. Also, we reduced the stell ar mass 
fraction by 30 per ce nt jLongjhetti fc Saraccol l2009l ) to ac- 
count for their use ofa lSalpetej(|l955h IMF. A dditionally, we 
include a result on Galaxy scales at z = from IConroy et al.l 
who determined stellar masses from the spectro- 
scopic SDSS value-added catalogue (|Blanton et al.l l2005h . 
along with halo mass estimates utilising the measured ve- 
locity dispersions of the satellites (we have scaled their halo 
mass, from M200 to M5 00, assuming an N FW profile with 
concentration given by iDuffv et al.l 20081 ') . We further as- 
sume that the stars are significantly more concentrated than 
the DM and that modifying the mass cut from -R200 to 
i?500 will have a negli gible effect on the stellar mass. We 
note that the result of Conrov et al.l (|2007l ) is in agreement 



ing to determine halo masses. 

The trend seen in our simulations, is a stellar frac- 
tion that decreases gently as M500 increases from 10^^ to 
10" Mp) . This is in accord w ith the semi-analytic model re- 
sults in iBalogh et al.l (|2008l ). Simulations with weak or no 
feedback contain stellar masses that are 2-3 times higher 
than observed, PrimC_WFB is an exception. This differ- 
ence is expected as the dominant cooling channel at the 
virial temperatures of haloes with mass 10" Mq is be- 
lieved to be from metal-line emission and hence 'ZC will 
have enhanced cooling rates over 'PrimC. The simulations 
with stronger feedback, ZC_SFB and ZC_WFB_AGN, are in 
better agreement with the observations. 

We can see that simulations which have relatively ef- 
ficient supernova feedback in the dense regions of haloes, 
or include an additional heating source in the form of an 
AGN, predict more realistic stellar fractions in present day 
Groups and Glusters with M500 ^ 10^^ M©. If we include 
the full mass range at low redshift, we can conclude that 
AGN feedback is necessary to prevent stellar mass building 
up in Gala xies, Groups and Gluste rs. This is in good agree- 
ment with [McCarthy et al.l |200^ , who found that model 
ZC_WFB_AGN predicts group K-band luminosities and gas 
fractions that are in excellent agreement with observations. 
However, in the following section we will see that at high red- 
shift the baryon-dominated inner halo of the schemes with 
weak or no feedback is required to match strong lensing re- 
sults, leading to an intriguing conflict which we will discuss 
further in Section [HI 
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Figure 3. Inner slope of the total mass density profile versus cen- 
tral baryon fraction. Only our results for Galaxy haloes at ^ = 2 
arc shown, t o compare with the obs ervational constraint on lens- 
ing galaxies jKoopmans et aTll2006l : the vertical size of the, lower, 
maroon box indicates the intrinsic la Gaussian spread across all 
redshifts while the horizontal size of the box is the 68 per cent 
confidence interval of the quoted stellar fractions, estimated by 
bootstrapping. Note that for the sample of early type galaxies 
the stellar fraction is essentially the baryon fraction). The differ- 
ent symbol sizes denote different mass ranges. Symbol type and 
colour are used to distinguish different simulations. The black, 
hatched region indicates the quartile spread of the DMONLY sim- 
ulation. As is clear from this figure, only simulations with high 
central baryon fractions reproduce the observed steep inner den- 
sity profiles of high redshift Galaxies. This is in contrast with the 
preferred simulation schemes at higher masses and lower redshift 



4.2 Observed inner profile slopes 

Gravitational lensing of light by an intervening galaxy en- 
ables us to probe the inner mass profile of the lens galaxy's 
halo. The slope of the inner mass p rofile for a sample of len s 
galaxies at 2 < 1 was presented bv lKoopmans et al.l (|2006t ). 
Surprisingly, the inner slope is strongly constrained to be 
close to isothermal (/3 = ~ —2) with no evidence for 

evolution. The region where the slope is measured (the Ein- 
stein radius is typically around 3h~^kpc) is comparable to 
the scales accessible in our highest-resolution simulations at 
z — 2, allowing us to compare the two results. 

Fig. [3] shows the inner slope of the total mass density 
profile versus central baryon fraction for our Galaxy haloes 
at jz = 2. When feedback is weak or absent, the inner slope 
is close to the isothermal value (/3 = —2); when feedback 
is strong, the slope is flatter and close to the DMONLY 
values (P ~ —1.4; shown as the black, hatched region, Sec- 
tionl5.3p. Comparing t hese results with the observations of 
iKoopmans et all (|2006l ). shown as the lower maroon box in 
the figure, only the simulations with weak or no feedback 
produce similar values to the o bservations. The iso t herma l 
profile at z = 2 is also seen in iRomano-Diaz et al.l (|2008h , 
who found that for their simulations the influx of subhaloes 
at late times, 2^1, acted to flatten the profile. 

This conclusion is apparently at odds with what can be 
drawn from the observed stellar fractions. On the one hand, 
strong feedback is required to keep cooling under control, 
such that the observed stellar mass fractions in groups and 
clusters at low redshift are reproduced. On the other hand, 
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Figure 4. Wc plot the mean density profile of a number of species, 
normalised by the average virial mass and radius, and scaled by 
to reduce the dynamic range. Here we can see the significant 
dynamical influence of the baryons in this particular object. This 
is a galaxy-sized object at z = 2 taken from ZC_WFB, averaged 
over 20 relaxed haloes, from the simulation with comoving length 
25 Mpc. The yellow triangles are the total matter profile; the 
black squares are the DM profile; the blue pentagons indicate the 
gas profile and the green circles denote the stellar profile. Error 
bars illustrate the 68 per cent confidence limits within each radial 
shell, estimated by bootstrap analysis. Vertical lines indicate the 
region where the inner slope is measured (0.025 < r/i?vir < 0.05). 
The arrow represents iJpo3, taken to be the resolution limit. The 
solid red curve is the best-fit NFW profile to the DM, over the 
outer region (r/iJvir > 0.05). The DM mass was divided by (1 — 
fh(T < iJvir)) such that the total mass of the DM component 
would equal Mvir- With the halo mass known the NFW profile 
has only one free parameter. The NFW curve has had this factor 
removed for this plot. The legend contains the mean virial mass 
and radius of the haloes, Mvir and i?vir respectively, the best- 
fit NFW concentration, Cvir, and inner profile slope, /3. Baryons 
strongly influence both the total and the DM density profiles 
within 0.1 i?vir- 



the feedback has to be weak (or absent), to generate the 
steep inner profiles observed in galaxies at low and interme- 
diate redshift (out to z ~ 1). 

One possible reason for this dichotomy, is that the sim- 
ulated galaxies are at higher reds hift than the l ensing ob- 
servation^. However, as shown bv lDehnenI (|2005l ). the den- 
sity profile of a colhsionless system does not steepen during 
mergers, so the lensing galaxies must have been isothermal 
at higher redshift, where it is more likely that significant gas 
condensation could have occurred. Secondly, we cannot rule 
out selection effects in the observations that may have a bias 
towards steeper density profiles. This could be important 
as isothermal profiles do exist in the ZC_WFB_AGN sim- 
ulation, albeit significantly fewer in number than ZC_WFB 
creates. A final possibility is that the simulations themselves 
do not yet accurately model the high-redshift growth of the 
brightest galaxies. One would therefore require a feedback 
prescription that was less effective at high redshift (allow- 
ing more gas to accumulate and steepen the central profile). 



^ We use the highest resolution simulation at 2 = 2 due to 
th e close match b e tween typical Einstein radius of galaxies 
in iKoopmans et al. I ||2006|'| . ^ 3 h ^ kpc, and our fitting range 
Ri2-4.5/i-lkpc. 



The feedback must then rapidly expel the gas before it can 
form stars, faster than the dynamical timescale of the halo, 
such that the DM retains an isothermal profile. 



5 HALO DENSITY PROFILES 

To first get a general idea of the effects of baryons on the 
DM halo, we show an example of the halo density profile and 
its components in Fig. |4l where we have plotted the total 
mass profile, as well as the individual gas, stellar and DM 
components. Here, the halo is an average over 20 relaxed 
Galaxy haloes from ZC_WFB at z — 2, and is shown in 
dimensionless form: (r /Rvi^)'^${r / R^i^), where <I'(r-/7?vir) = 

p(r/i?vir)i??ir/Mvir. 

As can be seen in the figure, the baryons are more cen- 
trally concentrated than the DM. The stellar component has 
a steeper profile than the DM at all radii, whereas the gas is 
steeper in the inner region (r/i?vir < 0.05) and slightly fiat- 
ter in the outer region {r/Rvh > 0.1). The inner region of the 
DM profile has been significantly affected by the baryons. 
To highlight thi s change, we also sh ow the best-fit Navarro, 
Frenk & White (|Navarro et al.lll997l : hereafter NFW) profile 
(solid, red curve), fit to data with 0.05 ^ r/R^ii ^ 1. The 
NFW profile, which is a good approximation to the equiva- 
lent DMONLY profiles, is of the form 

Pcrit 

where Sc is a characteristic density contrast and is a scale 
radius. We fit for only one parameter, Vs, as we use the 
previously-measured virial radius and mass for each halo to 
define 5^. 

The simulated DM profiles are significantly steeper in 
the inner regions than the NFW profile. This is expected: 
the high central baryon fraction in this run (which, as we 
can see, is mainly in stars) has pulled the DM in towards 
the centre. Measurement of the inner slope of the DM profile 
(/3 ~ —2) confirms that the inner profile is isothermal. It is 
therefore clear that the effects of baryonic physics (radiative 
cooling in this case) can have a profound impact on the inner 
DM density profile. 



(r/r-,)(l + r-/r,)2 



(5) 



5.1 Galaxies and groups at high redshift 

We present mean DM density profiles for the z = 2 data in 
Fig. [51 These consist of 32 haloes from the 25 h^^Mpc simu- 
lation at Dwarf Galaxy and Galaxy mass scales (top panels) 
and an additional 14 Group scale objects (bottom panels) 
drawn from the larger 100/i~^Mpc simulation. Note that 
the group scale objects do not satisfy our stringent criterion 
for the inner profile study: the black arrow, depicting the 
P03 resolution limit, is at r/R^ir ~ 0.04. 

We have grouped the simulations together according to 
the physics prescriptions that we are testing. The first set 
is shown in the left panels and corresponds to the introduc- 
tion of metal-line cooling and weak stellar feedback. The 
second group, in the right panels, tests the types of feed- 
back, namely stellar (both weak and strong) and AGN, all 
with the same cooling prescription. The DM density profiles 
have been divided by (1 — /b"'^) for comparison with the 
DMONLY profiles. 
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Figure 5. In order to quantify the response of the DM halo to the varying physics prescriptions, we plot the mean DM halo density 
profile of haloes in fixed mass bins, Mvir = 5 X lO" - 5 X lO^'^h'^ Mq, the mean mass of which corresponds to Galaxies (top panel) and 
^vir = 10^'^ — 5 X 10^^h~^ Mq (the latter limit is to include the largest object) which are Groups (bottom) at ^ = 2 (drawn from the 25 
Mpc simulations). Units are normalised by the virial mass and radius from the equivalent haloes in the DMONLY simulation. The 
DM density profiles for the models including baryons have been divided by (1 — Z^"'^) to facilitate comparison with DMONLY. The left 
panels show results for the runs with weak or no feedback and the right panels for different feedback schemes. The vertical lines denote 
r/Rvir = 0.025 and 0.05, the region where the inner profile slope is estimated. The vertical arrows denote the P03 resolution limit; in 
the case of the Groups, the inner profile slope cannot be reliably measured at z = 2 (the other haloes have a resolution limit within the 
innermost line). Error bars are bootstrap estimates of the 68 per cent confidence limits about the mean density within each bin. The 
legend contains the mean virial mass of the haloes, Mvir, the best-fit NFW concentration, Cvir, and inner profile slope, /3. 




Figure 6. As in Fig. [5]but now for Groups and Clusters at z = (drawn from the 100 /i~^Mpc simulations). The top panels use a mass 
range Mvir = 10^^ — 10^'^h~^ Mq and the bottom panel for the range Mvir = 10^'* — 6 X lO^^/i"^ Mq (the upper bound includes the 
most massive cluster in the simulation). 
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As expected, the largest differences occur in the in- 
ner regions, where the density is highest. In the left pan- 
els, showing results for runs in which feedback effects are 
weak or absent, there is a significant steepening of the 
density profile towards higher densities on smaller scales 
across all mass ranges. The curves then more-or-less con- 
verge with DMONLY at radii larger than r/Rviv = 0.2. 
Given that PrimC.WFB and PrimCNFB predict nearly 
identical DM profiles, but differ strongly from DMONLY, it 
is clear that cooling plays a crucial role in determining the 
magnitude of the back-reaction on the DM. Indeed, com- 
paring PrimC_WFB with ZC_WFB, we see that including 
metal-line cooling results in a steeper profile for Galaxies. 
This difference is reduced in more massive systems for which 
the virial temperatures exc eed the regime in wh ich metal- 
line cooling is efficient (e.g. IWiersma et al.ll2009l ). 

In the right panels of Fig. [S] we see that the stronger 
feedback schemes (ZC.SFB and ZC.WFB_AGN) have DM 
density profiles that are closer to those from DMONLY. 
These runs have smaller baryon fractions, thus the overall 
effect of the cooling has been reduced. 

5.2 Groups and clusters at the present day 

For Groups and Clusters at 2: = we repeat our previous 
investigation in Fig. [S] Again, the DM profiles in the simu- 
lations with baryons diverge from the DMONLY results in 
the inner regions, although the effect is not as dramatic as it 
was for the less massive systems shown in Fig. O This is be- 
cause the typical cooling times are longer in these systems. 
Furthermore, the difference between the runs with weak and 
strong steUar feedback (ZC.WFB and ZC_SFB; right panels) 
is smaller, because the strong stellar feedback is less effec- 
tive in the more massive Groups and Clusters than it was in 
the lower-mass systems at high redshift. The ZC_WFB_AGN 
model (with the lowest central baryon fractions) produces 
almost identical profiles to DMONLY in Clusters, but a 
slightly lower profile in Groups. 

Another interesting effect apparent in Fig. [6] for the 
Clusters at z = concerns the difference between the 
DM profiles of the no and weak feedback schemes (left 
panels). The model showing the highest central DM den- 
sity is PrimC-WFB, while PrimC_NFB and ZC_WFB are 
statistically indistinguishable. At first glance this is unex- 
pected, as the latter two models have higher central baryon 
fractions (see Fig. [Tl b o ttom panel) . A similar result was 
found by iPedrosa et all (|2009l ) for their galaxy simulations 
at z = 0. We summarise their explanation of the effect within 
the context of our simulations. When weak stellar feedback 
is included (going from PrimC_NFB to PrimC_WFB), a cer- 
tain amount of the gas will be expelled from satellite galaxies 
but not from the main (group or cluster) halo itself. As a 
result, the satellite haloes will be less bound and suffer more 
tidal disruption as their orbit decays due to dynamical fric- 
tion. The result is that less mass is transferred to the centre 
of the halo. When metal enrichment is added to the simu- 
lation (going from PrimC.WFB to ZC.WFB), the coohng 
time of the satellite gas becomes shorter (due to metal-line 
emission) and as a result, the satellite is able to hold on to 
more of its gas. This reduces the effect the feedback has on 
the evolution of the satellite halo itself (and thus more mass 
can be transferred to the centre). The most important con- 



sequence of these effects is where the angular momentum of 
the satellites gets transferred. In the case of PrimC_WFB, 
less of the angular momentum is transferred to the inner 
region than in PrimC_NFB, for example. As a result, the 
inner profile is denser in the former case, even though the 
overall baryonic mass is smaller in the centre of the halo. 



5.3 Inner profiles 

It is clear that the main driver of the change in the (inner) 
DM profile is the condensation of gas to smaller scales than 
the DM, thereby increasing both the local baryon fraction 
and the DM concentration in the now-deepened potential 
well. If we assume that DM is pulled inwards by the condens- 
ing baryons, then we should expect some relation between 
the DM density profile and the baryon fraction within the 
inner region. We test this in Fig. [7] by plotting the median 
inner power law slope of the DM profile, /3dm, as a function 
of the median fh{r/Rvir ^ 0.05) value. The quartile spread 
in /3 values for the DMONLY run is also illustrated, as the 
hatched region. We remind the reader that in the inner pro- 
file study we utilised a well-resolved subsample of the full 
halo catalogue, such that -Rpoa ^ 0.025-Rvir. 

For the z = Group and Cluster haloes (bottom panel), 
/b varies by more than a factor of 3 in the inner region of the 
haloes and yet the resultant inner slope stays within 20 per 
cent of the DMONLY value. There is a tentative steepen- 
ing of the inner profile with increasing baryon fraction, but 
the steepest profile found in the simulation (PrimC_WFB) 
does not have the largest central baryon fraction, as already 
discussed in Section [5.21 

Our z = 2 Galaxy sample (top panel in Fig. [7| does 
demonstrate a clear and significant trend of steepening inner 
DM density profile with increasing central baryon fraction. 
The lowest baryon fractions are found in the ZC_WFB_AGN 
simulation for which the DM haloes are indistinguish- 
able from those in DMONLY. In models PrimC_NFB and 
ZC_WFB the baryon-dominated central regions generate 
nearly isothermal (/3 — —2) inner DM density profiles. 



6 HALO CONCENTRATIONS 

For an NFW halo of a given mass, the DM density pro- 
file is specified entirely by one parameter, the concentra- 
tion. In this section, we measure and compare NFW con- 
centrations for the DM profiles of haloes in our simulations 
(Section [GTIJ. We then consider the total mass profile of the 
halo by utilising simple, non-parametric measures of the con- 
centration, based on ratios of density contrasts at different 
scales (Section I6.2|l and the velocity profile of the halo, in 
particular the maximum velocity, in Section [6.31 All haloes 
with a P03 convergence radius, iipos < 0.05-Rvir, are con- 
sidered in this work on halo concentrations. This sets an 
effective mass limit of Mvir > 8 x 10^" Ii^^Mq at z = 2 and 
Afvir > 5 X 10^^ /i"^M0 at z = 0. The sample of haloes in- 
cludes over 500 objects at z = and ~ 20 (~ 150) systems 
at z = 2 in the simulation volumes with comoving length 
100 (25) ft-^Mpc. 
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Figure 7. The median inner (0.025 ^ f/Kvir ^ 0.05) power 
law slope of the DM density profile as a function of the median 
baryon fraction within r = 0.05iJvir! for different simulations at 
z = 2 (top) and 2 = (bottom). The symbols (and colours) rep- 
resent different simulations while the symbol size indicates the 
mass range. The vertical error bars illustrate the quartile scatter. 
The black hatched region represents the quartile spread of the 
DMONLY simulation. Generically, higher central baryon concen- 
trations yield steeper inner DM density profiles. 



6.1 DM profile: NFW concentrations 

A well established result from A^-body simula- 
tions (|Bullock et al.llioOll 'l is that the NFW concentration of 
the DM halo is anti-correlated with its m a,ss (fo r the latest 



results see iNeto et'Z' '200'^; iDuffv et al.) l2008l : iGao et all 



|2008| and iMaccio ct al . 20^) and takes on a power-law 
form 



(6) 



where Bvir is close to —0.1 when fit to data over nearly 
five orders of magnitude in mass. However, it is not clear 
how much this trend, which is primarily driven by the for- 
mation time of the halo, is modified by the presence of 
bary ons. Observations of X-ray luminous group s and clus- 
ters iBuote et all l2007l : ISchmidt fc AUenl I2OO7I ') suggest a 
steeper dependence of concentration on mass t han the DM 
only s imulations predict, as was pointed out bv lDuffv et al.l 
(|2008l ) . This is primarily due to observed groups having ~ 30 
per cent higher concentrations than the simulated objects 
(the concentrations of clusters were in good agreement with 
the simulations if a subsample of dynamically-relaxed haloes 



Figure 8. We plot the NFW DM halo concentrations from the 
baryon simulations, normalised by the best fit concentration-mass 
relation from DMONLY, as a function of halo virial mass at 2 = 2 
(0) in the top (bottom) panel. Values greater than unity indi- 
cate that the DM halo has contracted under the influence of the 
baryons. The points represent the median concentration within 
each mass bin. The vertical error bars are the 68 per cent con- 
fidence intervals estimated by bootstrapping the haloes within 
each mass bin. The horizontal error bars indicate the mass range 
of each bin. Note that the gap in the mass coverage at 2 = 2 is 
where the two simulation volumes meet. In the absence of strong 
feedback and when metal-line cooling is included, baryons sub- 
stantially increase the NFW DM concentrations of Galaxies, but 
the effect on Groups and Clusters is much smaller. AGN and 
strong supernova feedback actually reduce the concentrations of 
Groups. 



was used in the comparison). It is therefore important to 
check whether the inclusion of baryons can bring theory and 
observations into agreement. 

We will present concentrations relative to the equiv- 
alent values for the DMONLY model. Note that our sim- 
ulations assume the WMAP3 cosmology, which has an 8 
per cent lower value of a s than the WMAP 5 value (0.74 
versus 0.796) assumed bv iDuffv etHI (|2008l ') and leads to 
somewhat smaller concentrations at fixed masfl As de- 
scribed in Section 12.21 we fit density profiles over the range 
0.05 ^ r/Rvir ^ 1 We first assess the goodness-of-fit of the 

* For reference, including all resolved haloes, our best-fit power 
law relations for DMONLY are [Avir, Bvir] = [7.8 ± 0.6, -0.10 ± 



0.03] at 2 = 0, and [A^i^, 



[3.7 ±0.2, 0.01 ±0.03] at 2 = 2. 
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NFW function when baryons are included. To do this, we 
compute the usual quantity 



2 



iVbin 



E 

j=i 



(logio Psim.i - logio pNFW,i) 



(7) 



where A''bins is the number of bins in our profile and psim and 
Pnfw are the densities from the simulation and the best- 
fit NFW profile respectively. For DMONLY we find that 
o"flt ~ 0.02. For the runs with baryons the goodness-of-fit is 
similar (typically within f per cent) for Group and Cluster 
haloes. For smaller objects the difference is more pronounced 
for the simulations with high central baryon fractions, for 
which cTflt increases by around a factor 2. 

Shown in Fig. [S] are the NFW concentrations of the 
DM haloes in the runs with baryons, relative to the corre- 
sponding values from DMONLY, as a function of halo virial 
mass. For the Group and Gluster haloes at z = (bot- 
tom row), the only simulation to show substantial (> fO 
per cent) deviations from the DMONLY values, is the sim- 
ulation that includes AGN feedback, ZC_WFB_AGN. For 
this model the NFW concentrations of Mvir ~ 10^^ fe'^M© 
haloes are about 15 per cent lower. For the run with strong 
stellar feedback, ZC_SFB, the decrease is about 10 per cent. 
In both cases the expulsion of gas has ca used the D M to 
expand relative to the DMONLY case fe.g. lHillsl[l980h . 

The effect of baryons on the NFW concentrations of 
z = 2 Dwarf Galaxy, Galaxy and Group haloes is shown in 
the top row of Fig.[8l As before, the differences are more dra- 
matic for these lower mass, higher redshift objects. In the 
runs without strong feedback the concentration increases, as 
expected. The increase is typically 10-20 per cent for Groups, 
but can be as large as 50 per cent (when supernova feed- 
back is absent) for Dwarf Galaxies. Th is dramatic increase 
is sim ilar in magnitude to that found bv lRomano-Di'az et al.l 
120091 ') for the concentration of a Mvir ~ 10^^ h'^M^ halo. 
As was the case at z = 0, in runs with effective feedback 
the presence of baryons makes little difference to the con- 
centration of the DM; the maximum effect being a decrease 
of ~ 15 per cent for Galaxies in ZC_WFB_AGN. 

We checked how the results change w hen only relaxed 
haloes are selected (as defin ed in Duffv et al..,2008i ). As was 
found in previous work (e.g. iDuffv et al. 20081 ') . the average 
concentration of the haloes increases, but the power-law re- 
lation between concentration and mass remains the same 
within the errors. 

We have also checked explicitly that the concentration 
of a halo increases with its central baryon fraction. This is 
indeed the case for all simulations and mass scales, except for 
the Group and Gluster haloes at z = 0, in the PrimC_WFB 
model. As discussed in the previous section, this run has 
haloes on these mass scales with anomalously high central 
DM densities (as compared with the run with primordial 
coolin g and no feedback, PrimC_NFB). 

In Duffv et al.l (|2008l ) we demonstrated that the inferred 
NFW concentrations of groups observed in X-rays at z = 
are 30 per cent more concentrated than predicted from 
DM only simulations. It was suggested that the inclusion 
of baryons would alleviate this discrepancy through the con- 
traction of the DM halo. As is clear from Fig.[8l however, the 
largest increase of any physics scheme is still less than 10 per 
cent. Moreover, for those schemes which reproduce the ob- 
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Figure 9. Here we plot the median ratio of the sphcrieal over- 
density radii K500 and K2500 ^ ^ function of halo mass, M500 
in this case, at z = 2 (top) and (bottom). This ratio is a use- 
ful non-parametric measure of the concentration of a halo. The 
vertical error bars are 68% confidence limits about the medians, 
estimated by bootstrapping the samples within each bin, and the 
horizontal values represent the mass range of the bin. Note that 
higher values of the ratio indicate lower concentrations. All sim- 
ulations show a positive trend with mass; this is indicative of the 
decreasing amount of baryons able to cool and condense as the 
virial temperature of the halo increases. The offset in the nor- 
malisation between simulation schemes is larger than the scatter, 
allowing the use of this ratio as a nonparametric concentration 
proxy. 



served z — stellar fractions (shown in Fig. [2]), the strength 
of feedback is such that we actually reduce the concentration 
of the Groups relative to the DMONLY simulation. As an 
additional issue, the trend of the concentration ratio with 
mass is positive in the strong feedback schemes which may 
further increase the disagreement with observations. Clearly, 
the inclusion of baryons does not resolve the problem. It is 
thus important to check whether observational biases or se- 
lection effects can account for the mismatch between theory 
and observation. 



6.2 Total mass profile: Non-parametric 
concentrations 

Although the NFW profile is a reasonable approximation 
to the DM profile at large radii {r/Rvh > 0.05, as shown 
in Fig. 31, the measurement of its concentration parameter 
requires the shape of the profile to be constrained over a 
range of scales, with accurate removal of the baryonic mass 
profile. A simpler, and thus more easily achievable method 
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from an observational point of view, is to consider the entire 
halo, DM and baryon components together, and measure the 
mass/radius ratio of a halo at two different spatial scales. 

In Fig. [9] we plot the ratio of two radii that are com- 
monly used by observers (e.g. in observations of X-ray 
groups and clusters), -R5oo/^250cEI, as a function of halo 
mass. All runs demonstrate a positive, albeit weak, depen- 
dence on mass with significant run-to-run variations in the 
normalisation. (Note that runs with higher central baryon 
fractions will typically have lower R500 / R2500 ratios, because 
the value of i?2500 grows as a result of the increased central 
mass.) The deviation from DMONLY is at the sub-25 per 
cent level and much smaller if the feedback is strong. The 
differences between the models are qualitatively similar to 
those for the NFW concentrations discussed in the previous 
section. Note, however, that contrary to the NFW concen- 
trations, the non-parametric total matter concentrations are 
never significantly reduced relative to the DMONLY case. 



6.3 Concentration measures: circular velocity 

A key, and relatively easily observed measure for the 
structure of a halo is the circular velocity, Vc{r) — 
{GM{< r)/rY^'^ , in particular its maximal value, Vmax. Like 
^?50o/fe500 , it can be more robustly determined than the 
NFW concentration. The creation of realistic velocity pro- 
files for galaxies has, however, been a long standing issue 
in simulations within the CDM paradigm. For A'^-body only 
simulations, the maximum velocity is well approximated by 
the ana lytic solution to the NFW profile fmax = vdr ~ 



am 

2.17rs) (|Cole fc Lacevlll996l : [Navarro et al.lll996f ). with the 
final result that Wmax ~ Kir = {GM^ir / RvirY^'^ ■ Observa- 
tions find that the maximum velocity in the disk is sim- 
ilar to the virial velocity (e.g. iDutton et al.l l2005l l . Typi- 
cally, however, when simulations include baryons, either ex- 
plicitly or through adiabatic contraction models, the haloes 
will have a maximum velocity t hat is a factor of two higher 
than the halo virial velocity (e.g. Navarro fc Steinmet j|200ol : 



iDutton et aIll2007l : [Pedrosa et al.ll2010l ). This increase in the 
velocity ratio is a consequence of the contraction of the halo 
in the presence of a significant central baryon component, 
with a large velocity ratio indicative of a high NFW concen- 
tration. 

Recently, IPedrosa et"aLl (|20ld ) investigated the circular 
velocity profile in a resimulation of a single high resolution 
galaxy for various physics implementations. They found that 
there was a positive correlation between the ratio 
and the 'sharpness' of the DM density profil^f], indicative 
of the contraction of the halo in the presence of baryons. 
Furthermore, effective feedback was necessary to obtain a 
realistic, i.e. low, velocity ratio in their galaxy. In addition 
to galaxies we have investigated a larger mass range, groups 
and clusters. We also lend statistical weight to conclusions 
concerning the importance of the various physics implemen- 
tations by examining a large sample of haloes. A more ex- 
tensive investigation of the rotation properties of the OWLS 



galaxies at z = 2, in addition to Wmax, is being performed in 
a separate study (Sales et al., in prep). 

In Fig. [To] we show the maximum circular velocity, in 
units of the virial velocity, as a function of halo mass at 2 = 2 
(0) in the top (bottom) panel for various implementations 
of the cooling and feedback prescriptions. The most striking 
result is the good agreement between the strong feedback 
schemes (and DMONLY) with Umax ~ Vvi-c at all masses 
and redshifts, and the significant offset for the other physics 
implementations. This divergence is reduced with increasing 
mass at all redshifts due to the strong anti-correlation of the 
velocity ratio with halo mass in the weak/no feedback runs. 
This is indicative of the reduction in gas cooling efflciency as 
the halo mass, and hence the virial temperature, increases. 

At 2 = 2 (top panel) there is a dramatic divergence in 
the velocity profile below 10^^ /i~^M0 between the runs with 
primordial cooling, due to the supernova feedback. When 
metal-line cooling is included, ZC_WFB, the divergence with 
the no feedback model is reduced; the overall effect of the 
enhanced cooling is therefore to counteract the gas removal 
efforts of the supernovae. 

At z = Q (bottom panel of Fig. llOf) we once again find 
the counterintuitive result that the model with weak feed- 
back, PrimC_WFB, is more concentrated in clusters (but not 
in groups) than the model without feedback PrimC_NFB. 
We also see the familiar effect of the metal-line cooling 
overcoming the weak supernova feedback when comparing 
PrimCNFB and ZC.WFB models. 

It is clear that strong feedback is necessary if one wishes 
to limit the effect of the baryons in increasing the maximum 
circular velocity of the halo above the virial velocity. 



7 ADIABATIC CONTRACTION 

Having demonstrated that baryons can significantly infiu- 
ence the DM density profile, we will now assess the degree 
to which this modification can be modelled as a secular adi- 
ab atic contraction of the D M halo. We will test the models 
of iBlumenthal et al.1 l| 19861 '). henceforth B86, and G04, and 
will make use of the publicly-available code CONTRA (G04). 

The B86 model for adiabatic contraction assumes that 
the DM halo is spherically-symmetric and that the particles 
are on circular orbits so that we can think of the halo as a 
series of shells that can contract but do not cross. Assuming 
that the baryons initially trace the DM halo density pro- 
file and then fall slowly (i.e., such that the mass internal to 
radius r changes slowly compared to the orbital period at 
r) towards the centre as they cool, we can compute the re- 
sponse of the dark matter shells. In that case the dark mat- 
ter particles conserve their angular momentum and hence 
rvc oc [rM(r)]^^^ is conserved, where M(r) is the total mass 
internal to r. We thus have 



Mdm,i(n)n 

1 J'univ 

^ Jh 



[Mdm,f(rf) + Mb,f(r-f)] n, 



(8) 



^ Typical values for these radii in terms of -Rvir are 0.5 (0.6) 
and 0.2 (0.3) for -R500 a^nd -R2500 at z = (2) respectively. The 
values will change by 10 per cent dependent on the simulation 
physics due to the baryonic back-reaction. 



^ The sharpness of the DM halo is characterised by the rate at 
which the density profile becomes shallower towards the halo cen- 
tre. T his is explicitly modelled in the Einasto profile llEinastol 
Il965h as an additional parameter to the NFW in the form of a 
rolling power law. 
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Figure 10. We plot the maximum circular velocity of the halo, 
relative to the circular velocity at the virial radius, Rvin as a 
function of virial mass at 2 = 2 (0) in the top (bottom) panel. 
This measure is sensitive to the degree of contraction of the halo 
in the presence of baryons. Each point corresponds to the me- 
dian value for haloes within a given mass bin. Vertical error bars 
are bootstrap estimates of the 68% confidence interval for each 
bin and the horizontal bars are the bin widths. A number of fa- 
miliar halo responses is easily seen in this figure, with galaxies 
more sensitive to the presence of the baryons than higher mass 
systems. The very high maximum velocity in PrimC_NFB indi- 
cates the importance of feedback in forming observed galaxies, 
whose maxim um circular velocity in the disk is similar to the 
virial velocity jPutton et al .120051 ). Furthermore, the sharp diver- 
gence between PrimC_NFB and PrimC_WFB at a critical mass 
of 10^^ /i~^M0 indicates the mass range at which the weak feed- 
back model becomes ineffective. Finally, the metal-line cooling 
has the overall effect of counteracting the supernova feedback in 
Galaxies and Groups at z = 2, seen by the agreement between 
PrimC^FB and ZC.WFB. The anti-correlation seen at 2: = 0, 
for the simulations PrimC-NFB and ZC.WFB, is indicative of 
the rising halo virial temperature restricting the efficiency of gas 
cooling. The strong feedback schemes are in close agreement to 
the DMONLY simulation predictions at all masses and redshifts, 
and are therefore necessary to recover observed values. 



where Mdm,i('"i) is the initial DM profile, Afb, £(?"£) is the final 
baryon profile, and rf is the final radius of the DM shell that 
was initially (i.e. before the baryons contracted) at ri. Note 
that we made use of the equality A^dm.fXn) = A/dm,i(ri), 
which holds because DM shells do not cross. Assuming 
that we know both Mh,!{ri) and A/dm,i(n), we can solve 
for n (and thus Afdm.iCn) = MdmC^f)) as a function of 
rf. In the following we will make the standard assumption 
that Afdm,i(ri) is given by the density profile (reduced by 
(1 - .fb"'")) of the corresponding DMONLY halo. 

G04 found that this method did not provide a satis- 
factory description of their numerical simulations. Because 
the DM particles are typically not on circular orbits, they 
suggested replacing Af(r)r by M{f)r, where f is the orbit- 
averaged radius, which they found can be approximated as 



i?vir^(r/J?vir 



(10) 



where A — 0.8 5 and w = 0.8 (c. f. the B 86 model which has 
A = w = 1). iGustafsson et all (|2006l ) extended this work 
and showed that the values of A and w that best describe 
the difference between simulations with and without baryons 
depend on halo mass and the baryonic physics implementa- 
tion. Furthermore, they showed that the best fitting values 
of A and w generally differ substantially from the values 
that provide a good fit to the actual mean radius f of the 
DM particle orbits. This suggests that while the introduc- 
tion of the two free parameters A and w enables better fits, 
the underlying model does no t capture the relev ant physics. 
Here, we extend the work by IGustafsson et all (j2006, ) to a 
larger range of mass scales and gas physics models. 



7.1 Best-fit G04 parameters 

We first compute the best-fitting parameters A and w for a 
subset of our simulated haloes with baryons to see if our 
simulations prefer a specific combination and, if so, how 
this compares with those found in previous work. We fo- 
cus on 3 models: the runs with no feedback (PrimC_NFB), 
weak stellar feedback (ZC_WFB) and strong stellar feedback 
(ZC_SFB). For each simulation we matched the haloes to the 
corresponding objects in DMONLY, producing an average 
density profile at each redshift. At z = 2 we consider haloes 
with Afvir = [5 - 50] X 10" /i"^Mq, while at 2; = we study 
the range Afvir = [3 - 60] x lO" h'^MQ. 

In Fig. [11] we show the distribution of aat values in 
the A — w plane, for each of the three simulation models at 
z = 2 (top panels) and z = (bottom panels). We calculate 
(jfit using equation (O, replacing pnfw with the appropri- 
ate adiabatically contracted density profile, over the range 

0.025 «; r/j?vir «; 1. 

The best-fit values of A and w depend strongly on both 
halo mass and on the baryonic physics. There is significant 
degeneracy between the parameters A and w, higher values 
of A can be compensated by higher values of w. The pa- 
rameters suggested by Go4ll and the original model by B86 
work best for simulation ZC_WFB, but even here they differ 



^ Note that we have rescalcd their value of A to our definition of 
the virial radius. 
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Figure 11. The goodness of fit, cflt, of the G04 model for adiabatic contraction as a function of the parameters A and w, in logarithmic 
units. The best-fit parameter combination is indicated by the red cross. The top and bottom rows show results at z = 2 and z = 
respectively. From left-to-right we show the results for the PrimCJNFFB, ZC_WFB and ZC_SFB simulations respectively. For each 
simulation the results are averaged over haloes matched to the DMONLY haloes with virial masses in the range 5 — 50 X 10^^ (3 — 60 X 10^'^) 
M0 for the 2 = 2 (z = 0) samples. The parameter values corresponding to the models of B86 (A = w = I, hence top right corner) 
and G04 are shown as asterisk and open diamonds, respectively. The two parameters are substantially degenerate. The best-fit values 
depend on halo mass, rcdshift, and on the implemented baryonic physics. 



substantially from ou r best-fit values. These findings are in 
good agreement with iGustafeson et all l|2006l ). 



7.2 Predicted DM density profiles 

The DM density profiles predicted by the adiabatic con- 
traction models are shown in Fig. 1121 When the parameter 
values suggested by G04 are used, the G04 model predicts 
similar profiles to B86, which do not describe the contracted 
DM profiles well. The models typically underestimate the 
DM density for r > 10~^i?vir and more so for the simula- 
tions with stronger feedback. If, on the other hand, we use 
the best-fit values of A and w for each simulation and halo 
sample then the predictions of the G04 method agree much 
better with the simulated profiles, but even in that case one 
would obtain a closer match to the actual density profile by 
neglecting adiabatic contraction for r > O.lfJvir- 

It is perhaps not too surprising that the models for adi- 
abatic contraction do not describe the simulations well. The 
assumption that the baryons initially trace the DM halo 
profiles is clearly violated in hierarchical models. Haloes are 
built by mergers of smaller progenitors, and cooling and 
feedback have already redistributed the baryons in these 
objects. Rather than contracting slowly as the gas cools, 
a large fraction of the ba ryons simply fall in cold (|Kav et al.l 
I2OOOI : iKeres et al.l [ioosi '). Moreover, in the stronger feed- 
back models a substantial fraction of the baryons is ejected. 
iTissera et all (|2009l ) recently demonstrated that the con- 
traction is manifestly not adiabatic. They find that the 
pseudo-phase space density relation is strongly modified 



when baryons were added to high resolution DM only sim- 
ulations. 

Models for adiabatic contraction are required when full 
hydrodynamic simulations are not available. Unfortunately, 
it is not possible to predict what values of A and w to use for 
the G04 model without a much better understanding of the 
baryonic physics. Moreover, even if the physics were known, 
we would need to simulate many haloes because the best- 
fit values of the parameters depend on both halo mass and 
redshift. 

Our results suggest that it is better to ignore adiabatic 
contraction for r > 0.1 Rvir, but that the use of adiabatic 
contraction models such as those by B86 and G04 does typ- 
ically represent an improvement at smaller radii, provided 
the feedback is moderate or weak. 



8 CONCLUSIONS 

Our main aim in this work was to investigate the response 
of the DM halo to the presence of baryons at a variety of 
masses and redshifts. We utilised a series of high-resolution 
simulations within a cosmological volume, with a number of 
different prescriptions for the sub-grid physics, to probe the 
effect of baryons on the DM distribution of haloes. Our re- 
sults centred on galaxy-scale haloes at z = 2 and groups and 
clusters at 2 = 0. We were particularly interested to discern 
the effect of the baryons when going from the situation in 
which radiative cooling dominates (leading to a high cen- 
tral concentration of baryons in the form of stars and cold 
gas) to one where feedback dominates (reducing the central 
galaxy mass and expelling gas from the halo). As we showed 
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Figure 12. We test adiabatic contraction models by comparing the predicted mean DM density profiles, in average virial units of the 
matched DMONLY haloes, for the simulations PrimC_NFB (left), ZC.WFB (middle) and ZCSFB (right) respectively. Also shown are 
the mean DMONLY profiles (which have been reduced by (1 — /^"'^) for comparison). Overlaid are the predictions from the adiabatic 
contraction models of B86 (triple-dot-dashed), G04 using their default parameter values (scaled to our definition of the virial radius) 
(dashed), and G04 with our best-fit parameter values (solid) as determined separately for each model, mass range, and redshift in Fig, 1111 
along with the goodness-of-fit measure (Tfit, Additionally, in the top legend we give the mean virial mass of the haloes, Mvir, the best-fit 
NFW concentration, Cvir, and inner profile slope, /3. The top (bottom) row shows results for haloes at 2 = 2 (^ = 0), No one parameter 
combination can reproduce the range of DM haloes, indicating that the slowly cooling gas picture, on which adiabatic contraction models 
are based, is not sufficient to model the behaviour of the DM in a live simulation. In fact, ignoring adiabatic contraction altogether gives 
the best results for r > 0,1 i?vir- 



in Fig, [T] when comparing central and global baryon frac- 
tions, our simulations with strong feedback (ZC_WFB_AGN 
and ZC_SFB) reduce the baryon fraction /b, in comparison 
with PrimC_NFB (no feedback) , by factors of 2-3 for Galaxy 
haloes at 2 = 2. In Group and Glusters the AGN can remove 
nearly half of the baryons from the inner region of the halo, 
at 2: = 0, 

By comparing with observed stellar fractions in low 
redshift Groups and Clusters of galaxies (Fig. [2]) we found 
that the simulations with a high baryon fraction (Fig, [l]) 
also predict stellar fractions significantly larger than ob- 
served. The simulation with inefficient gas cooling and stel- 
lar feedback, PrimC_WFB, and the strong feedback mod- 
els, ZC_WFB_AGN and ZCSFB, are broadly in agreement 
with the observed stellar fractions in 2: = objects of mass 
Mvir ~ 10^* /i~^M0, However, observed maximum star for- 
mation efficiencies of order 10% - 20% are only reproduced 
by the inclusion of AGN feedback. 

However, these same strong feedback simulations are in 
disagreement with the constraints inferred from combined 
gravitational lensing and stellar dynamics analyses of the 
inner total mass density profile of massive, early-type galax- 



ies (Fig, O , In this case the efficient feedback prevents the 
steepening of the density profile that is necessary to repro- 
duce the observed isothermal profiles. Instead, the simula- 
tions with high baryon fractions are in closer agreement, A 
more detailed comparison between the observations and sim- 
ulations is warranted, especially with regards to the biasing 
of lensing observations to steeper density profiles. 

An enhanced baryon fraction in the inner halo is ex- 
pected to contract the DM distribution, as a response to the 
deeper potential well of the system. This effect was clearly 
seen, from Dwarf Galaxy to Cluster scales and at both low 
(2 = 0) and high (2 — 2) redshifts, and is summarised 
in Fig, [T] where haloes with larger central baryon frac- 
tions develop steeper central profiles (especially the Galaxy 
haloes at high redshift). To quantify the contraction effect, 
we fit NFW profiles and compared them to a simulation 
with no baryons. Variations in the concentration are typ- 
ically around 20 per cent or less, although the concentra- 
tion of high-redshift dwarf galaxy haloes can be as much 
as 50 per cent higher when feedback effects are ignored. 
Strong feedback produces a mild decrease in the concentra- 
tion of a halo through the removal of a significant amount 
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of gas from the halo, in a similar manner to what was found 
bv lKovama et all (j2008h . 

We also investigated non-parametric measures of con- 
centration. We found that the ratio R^oo/ R2500 and the max- 
imum circular velocity are both useful indicators of the de- 
gree of contraction of the total mass profile. Efficient feed- 
back is required to redistribute the mass in the halo such 
that the maximum circular velocity is similar to the virial 
velocity, as is found observationally for disk galaxies. 

An interesting counterexample to the rule of increas- 
ing baryonic condensation leading to more concentrated DM 
haloes was witnessed. In our low redshift groups and clus- 
ters, the dark matter is denser in the simulation with weak 
stellar feedback (PrimC_WFB) than in the simulation with 
no feedback (PrimC_NFB), even though the central baryon 
fractions are lower in the latt er case. A similar result was 
found bv lPedrosa et al] (|2009l ). who analysed simulations of 
galaxies with and without feedback. They argued that the 
contraction of the DM halo is slowed by the infall of satellites 
(due to dynamical friction). If the satellites are less-bound, 
they may be tidally disrupted before they sink to the centre 
and their effect on halting the contraction is thus reduced. 

We compare d our results to the adiabatic contraction 
models of ^Blumcnthal et al.l (| 19861 ) and the revised model of 
iGnedin et all (|2004l ). We found that adiabatic contraction 
only improves the fits in the inner parts of the halo. Outside 
0.1-Rvir one would do better not to use any prescription at 
all. Within 0.1 i?vir the former model is unable to reproduce 
the back-reaction on the DM. The latter model can do a rea- 
sonably good job, but only if its two parameters are allowed 
to vary with the baryonic physics, halo mass, and redshif t, 
in agreement with the findings of iGustafsson et al.l (|2006h . 
This is of particular importance when feedback effects are 
strong, as required to reproduce the cooled baryon fractions 
in low redshift groups and clusters. This therefore removes 
the predictive power of the model and we caution its use, 
particularly if a detailed prediction for the structure of a 
DM halo is required. 

If one wishes to match the stellar fractions at low red- 
shift, then models with strong feedback are required to 
suppress star formation. The total mass profiles in such 
simulations are very similar to the DMONLY case, with 
iimax ~ Vvir, as obscrved. Intriguingly, the NFW concen- 
tration in these systems is actually reduced relative to the 
same halo in the DMONLY run. The finding reported in 
I Duffy et al.l (|2008l ') that observed groups appear to be signif- 
icantly more concentrated than simulated haloes, is there- 
fore still unexplained. The hypothesis that the inclusion of 
baryons would resolve the discrepancy between theory and 
observation has been shown to be wrong. Worse, the dis- 
agreement actually grows larger if one utilises strong feed- 
back physics schemes that can reproduce the observed stellar 
fractions in these systems. 
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APPENDIX A: RESOLUTION TESTS 

To assess the effects of resolution on our results, we compare 
the 512^ DMONLY and ZC_WFB simulations with lower 
resolution versions (256'^ and 128^, with softening lengths 
increased by factors of 2 and 4 respectively). We initially se- 
lect haloes with 10^ DM particles from the lowest resolution 
simulation. Matching haloes in the higher resolution simu- 
lations are found by first identifying high-resolution haloes 
that lie within the virial radius, i?po?, of the low resolution 
halo. We consider a halo to match if the candidate object 
fulfils the following criteria 

^^^^0.2; ^--5- ^0.2. (Al) 

This final selection removes ~ 10 per cent of the haloes from 
the sample, these objects are in the process of merging and 
hence would create artifacts in the averaged density profiles. 
We therefore perform resolution tests on 57 (27) matched 
haloes from ZC_WFB at z = (2) and 51 (27) haloes from 
DMONLY at 2 = (2). Note that this matching is different 
to the scheme employed in the rest of this work which was 
based on linking identical DM particles between different 
physics simulations. 

Fig. I A 1 1 compares average total mass density profiles for 
a group and a cluster-scale halo at 2; = in the three runs. 
In the DMONLY case, resolution effects cause the density to 
be underestimated, whereas the density is overestimated for 
ZC_WFB. At z = 2, the results from (lower mass) dwarf and 
galaxy-scale haloes are similar, although the density is also 
underestimated in the under-resolved regions in ZC_WFB, 

(Fig..m , „ 

iPower et all (|2003l ') defined a convergence radius, such 

that the mean two-body relaxation time of the particles 
within that radius is of the order a Hubble time. The equiv- 
alent radius is shown for each of our simulations in the two 
figures, as a solid vertical arrow. Although this radius was 
originally applied to coUisionless Ai'-body simulations (like 
DMONLY), we can see that it also provides a satisfactory 
indication of where the density profile becomes numerically 
resolved if baryons are included. This is especially visible in 
the lower sub-panels, where the fractional deviations of the 
two lower-resolution profiles from our standard-resolution 
(512^) results are shown. 
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Figure Al. Comparison of moan DM density profiles (and their residuals) for Group and Cluster haloes (top and bottom panels 
respectively), drawn from the 100/i~^Mpc box at z = 0, for runs with different resolutions. Results from DMONLY are shown in the 
left column, while results from ZC_WFB are shown in the right column. The colours indicate the resolution, from 128^ (red), to 256^ 
(blue) to 512^ (black). Most softening lengths are sufficiently small to fall outside of the plotted area but, where visible, softening scales 
are indicated with vertical lines with colour and line-style indicating simulation resolutions, described in the legend. The arrows indicate 
the P03 convergence radius for each simulation, according to colour. The vertical lines denote 2.5 and 5 per cent of the virial radius, 
corresponding to the scale within which the inner density profile slope is measured. The legend contains the mean virial mass of the 
haloes, Mvir. 
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